var Yt Ct Bt Lt Lpt Kt Nt Wt Rt Qt Zt Gt Ot Mut;

parameters beta varphi delta psi F nu alpha theta phi xis xie;


load xpar;
 
beta       = xpar(1);
varphi     = xpar(2);
delta      = xpar(3);
psi        = xpar(4);
F          = xpar(5);
nu         = xpar(6);
alpha      = xpar(7);
theta      = xpar(8);
phi        = xpar(9);
xis        = xpar(10); 
xie        = xpar(11);

external_function(name = interpv,  nargs = 2, first_deriv_provided, second_deriv_provided); 
 

model;

1          = beta*Ct/Ct(+1)*(Rt(+1) + 1 - delta);
F*Wt       = (1 + xie)*beta*Ct/Ct(+1)*Qt(+1);
Qt         = Yt*Gt + beta*(1 - varphi)*Ct/Ct(+1)*Qt(+1);
Wt         = psi*Ct*Lt^nu; 
Ot         = (phi*((Rt/alpha)^alpha*(Wt/(1 - alpha))^(1 - alpha))^(1 - theta) + (1 - phi))^(1/(1 - theta)); 
Yt         = Ct + Kt - (1 - delta)*Kt(-1) + Bt; 
Wt         = (1 - alpha)*phi/Mut*((Rt/alpha)^alpha*(Wt/(1 - alpha))^(1 - alpha)/Ot)^(1 - theta)*Yt/Lpt; 
Rt         =       alpha*phi/Mut*((Rt/alpha)^alpha*(Wt/(1 - alpha))^(1 - alpha)/Ot)^(1 - theta)*Yt/Kt(-1);
1          = (1 - phi)/Mut*(1/Ot)^(1 - theta)*Yt/Bt;

Mut        = interpv(Nt(-1), 1)/(1 + xis);
Zt         = interpv(Nt(-1), 2);
Gt         = interpv(Nt(-1), 3)*(1 + xis);

Ot         = Zt/Mut; 

Lt         = Lpt + F*(Nt - (1 - varphi)*Nt(-1));

end;


load x1_ss;
load x2_ss;


initval;

Kt   = x1_ss(1); 
Nt   = x1_ss(2);

end;


endval;


Yt   = x2_ss(1);
Ct   = x2_ss(2);
Bt   = x2_ss(3);
Lt   = x2_ss(4);
Lpt  = x2_ss(5);
Kt   = x2_ss(6);
Nt   = x2_ss(7);
Wt   = x2_ss(8);
Rt   = x2_ss(9);
Qt   = x2_ss(10);
Zt   = x2_ss(11);
Mut  = x2_ss(12);
Gt   = x2_ss(13);  
Ot   = Zt/Mut; 

end;

resid;
%steady;

perfect_foresight_setup(periods = 250);
perfect_foresight_solver(tolf = 1e-9, tolx = 1e-9);
